!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
program main
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_amr
   use interface_block_destroy
   use interface_block_output
   use interface_model_initial
   use interface_load_state
   use interface_save_state
   use interface_memory_usage
   implicit none

   character(4) :: nt_char, cpu_char

   real(8) :: t_start, t_stop, t_wall
   real(8) :: max_speed_global

   type(block), pointer :: base, current

   integer(4) :: output_method, output_level, format, i_level
   integer(4) :: output_time_start, output_time_end
   character(128) :: datavis_path, var_name_list
   character(256) :: header
   namelist /datavis_options/                                                                                                      &
      output_method,                                                                                                               &
      output_level,                                                                                                                &
      output_time_start,                                                                                                           &
      output_time_end,                                                                                                             &
      datavis_path,                                                                                                                &
      var_name_list,                                                                                                               &
      format,                                                                                                                      &
      header

!----------------------------------------------------------------------------------------------------------------------------------|
!     mpi initialization
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_init(error)
   call mpi_comm_size(mpi_comm_world, cpu_number, error)
   call mpi_comm_rank(mpi_comm_world, cpu_id, error)
   cpu_id_p1 = cpu_id + 1

!----------------------------------------------------------------------------------------------------------------------------------|
   t_start = mpi_wtime()
   t_loop_e = t_start

!----------------------------------------------------------------------------------------------------------------------------------|
   if (cpu_id .eq. 0) then
      write(*, *) 
      write(*, '(a)') '|---------------|'
      write(*, '(a)') ' Program starts:'
      write(*, '(a)') '|---------------|'
      write(*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   t = 0.0d0
   t_hc = 0.0d0
   t_wall = 0.0d0

   call getcwd(directory)
   directory = trim(directory) // '/'

!----------------------------------------------------------------------------------------------------------------------------------|
   t_s = mpi_wtime()

   call read_parameters
   if (cpu_id .eq. 0) call system("mkdir -p " // trim(datamap_path))
   call mpi_barrier(mpi_comm_world, error)
   
!----------------------------------------------------------------------------------------------------------------------------------|
   open(unit_datamap_manage_module, file = 'data_manage_precision.include')
   write(unit_datamap_manage_module, '(a, i)') 'integer(4), parameter :: output_real8_flag = ', output_real8_flag
   close(unit_datamap_manage_module)
   open(unit_datamap_manage_module, file = 'data_manage_parameter.include')
   write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(parameters_filename)),                              &
      ') :: parameters_filename = ''', trim(parameters_filename), ''''
   write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(datamap_path)),                                     &
      ') :: datamap_path = ''', trim(datamap_path), ''''
   close(unit_datamap_manage_module)

   open (unit_datamap, file = trim(parameters_filename))
   read (unit_datamap, nml = datavis_options)
   close (unit_datamap)
   open(unit_datamap_manage_module, file = 'data_manage_var_define.include')
   if (output_real8_flag .eq. 4) then
      write(unit_datamap_manage_module, '(a)') '   typedef float VAR_TYPE;'
   elseif (output_real8_flag .eq. 8) then
      write(unit_datamap_manage_module, '(a)') '   typedef double VAR_TYPE;'
   endif
   write(unit_datamap_manage_module, '(a, i, a)') '   int output_real8_flag = ', output_real8_flag, ';'
   write(unit_datamap_manage_module, '(a, i, a)') '   int output_time_start = ', output_time_start, ';'
   write(unit_datamap_manage_module, '(a, i, a)') '   int output_time_end = ', output_time_end, ';'
   write(unit_datamap_manage_module, '(a, a, a)') '   char datavis_path[128] = "', trim(datavis_path), '";'
   write(unit_datamap_manage_module, '(a, a, a)') '   char var_name_list[128] = "', trim(var_name_list), '";'
   write(unit_datamap_manage_module, '(a, a, a)') '   char datamap_path[128] = "', trim(datamap_path), '";'
   close(unit_datamap_manage_module)

!----------------------------------------------------------------------------------------------------------------------------------|
   allocate (base)
   allocate (base%x(nx_total), base%y(ny_total), base%z(nz_total),                                                                 &
      base%child(ini_block_x, ini_block_y, ini_block_z), base%p_pos(n_dims), base%id(n_dims, n_levels))

   call model_initial(base)

   call mpi_allreduce(max_speed, max_speed_global, 1, mpi_real8, mpi_max, mpi_comm_world, error)
   max_speed = max_speed_global

   base%dt = cfl / (max_speed / base%dx + max_speed / base%dy + max_speed / base%dz)
   base%dt = min(cfl * minval((/ base%dx, base%dy, base%dz /)) ** 2 / (max_et + eps), base%dt)
   base%dt_hc = min(cfl * minval((/ base%dx, base%dy, base%dz /)) ** 2 / (thermal_amp * max_te_ro + eps), base%dt) + eps
   min_dt = base%dt
   min_dt_hc = base%dt_hc

   if (load_state_flag .eq. 1) call load_state(base)

   t_e = mpi_wtime ()
   time_other = time_other + t_e - t_s

!----------------------------------------------------------------------------------------------------------------------------------|
   t_s = mpi_wtime ()

   if (cpu_id .eq. 0) then
      open(unit_par, file = trim(datamap_path) // 'params.mapout')
      write(unit_par, '(2x, a)') trim(directory)
      write(unit_par, '(a24, t32, i5)') 'cpu_number:', cpu_number
      write(unit_par, '(a24, t32, i5)') 'n_dims:', n_dims
      write(unit_par, '(a24, t32, i5)') 'sys:', sys
      write(unit_par, '(a24, t32, i5)') 'ng:', ng
      write(unit_par, '(a24, t32, i5)') 'ini_block_x:', ini_block_x
      write(unit_par, '(a24, t32, i5)') 'ini_block_y:', ini_block_y
      write(unit_par, '(a24, t32, i5)') 'ini_block_z:', ini_block_z
      write(unit_par, '(a24, t32, i5)') 'nx_no_bnd:', nx_no_bnd
      write(unit_par, '(a24, t32, i5)') 'nx_block:', nx_block
      write(unit_par, '(a24, t32, i5)') 'nx_total:', nx_total
      write(unit_par, '(a24, t32, i5)') 'ny_no_bnd:', ny_no_bnd
      write(unit_par, '(a24, t32, i5)') 'ny_block:', ny_block
      write(unit_par, '(a24, t32, i5)') 'ny_total:', ny_total
      write(unit_par, '(a24, t32, i5)') 'nz_no_bnd:', nz_no_bnd
      write(unit_par, '(a24, t32, i5)') 'nz_block:', nz_block
      write(unit_par, '(a24, t32, i5)') 'nz_total:', nz_total
      write(unit_par, '(a24, t32, i5)') 'ratio:', ratio
      write(unit_par, '(a24, t32, i5)') 'n_levels:', n_levels
      write(unit_par, '(a24, t32, es16.6)') 'gm:', gm
      write(unit_par, '(a24, t32, es16.6)') 'rangex:', rangex
      write(unit_par, '(a24, t32, es16.6)') 'rangey:', rangey
      write(unit_par, '(a24, t32, es16.6)') 'rangez:', rangez
      write(unit_par, '(a24, t32, es16.6)') 'zero_posx:', zero_posx
      write(unit_par, '(a24, t32, es16.6)') 'zero_posy:', zero_posy
      write(unit_par, '(a24, t32, es16.6)') 'zero_posz:', zero_posz
      write(unit_par, '(a24, t32, es16.6)') 'cfl:', cfl
      close(unit_par)
   endif

   t_e = mpi_wtime()
   time_output = time_output + t_e - t_s

!----------------------------------------------------------------------------------------------------------------------------------|
   do while (.true.)

!----------------------------------------------------------------------------------------------------------------------------------|
!     remaining time
!----------------------------------------------------------------------------------------------------------------------------------|
      t_loop_e = mpi_wtime ()

      t_loop = t_loop_e - t_loop_s

!      if (ns .gt. 10 .and. cpu_id .eq. 0) then
!         if (t_loop .gt. t_loop_pre * 4.0d0) then
!            write(*, *)
!            write(*, *) 'Warning: expend too much time in one loop.'
!            write(*, *) 't_loop_pre: ', t_loop_pre
!            write(*, *) 't_loop: ', t_loop
!            write(*, *) 'Minimum number of blocks = ', minval (n_blocks)
!            write(*, *) 'Maximum number of blocks = ', maxval (n_blocks)
!            write(*, *) 'time_regrid one loop = ', time_regrid_one
!            write(*, *) 'time_integrate one loop = ', time_integrate_one
!            write(*, *) 'time_boundary one loop = ', time_boundary_one
!            write(*, *) 'time_output one loop = ', time_output_one
!            write(*, *) 'time_other one loop = ', time_other_one
!            write(*, *)
!         endif
!      endif

      t_loop_pre = t_loop
      t_remain = abs(t_end - t) / (min_dt / t_loop)
      t_loop_s = t_loop_e

      time_other = time_other + time_other_one
      time_output = time_output + time_output_one
      time_boundary = time_boundary + time_boundary_one
      time_regrid = time_regrid + time_regrid_one
      time_integrate = time_integrate + time_integrate_one
      time_regrid_one = 0.0d0
      time_integrate_one = 0.0d0
      time_boundary_one = 0.0d0
      time_output_one = 0.0d0
      time_other_one = 0.0d0

!      if (t .ge. nt * dt_output .or. ns .eq. 99) then
      if (t .ge. nt * dt_output) then

!         call memory_usage(base)
!      call check_info ('output')

         t_s = mpi_wtime()

         ab_density = 0
         ab_pressure = 0
         nt = nt + 1

         if (cpu_id .eq. 0) then
            write(*,"(a, i7, a, i6, a, es10.3, a, es10.3, a, es10.3)")                                                        &
               ' ns= ', ns, ', blocks= ', n_blocks_all, ', t= ', t, ', dt= ', min_dt
         endif

         write(nt_char, '(i4.4)') nt - 1
         write(cpu_char, '(i4.4)') cpu_id

         if (cpu_id .eq. 0) then
            open(unit_nt, file = trim(datamap_path) // 'nt.mapout', form = 'unformatted')
            write(unit_nt) nt
            close(unit_nt)
            open(unit_t, file = trim(datamap_path) // 'info.t' // nt_char // '.mapout', form = 'unformatted')
            write(unit_t) t
            write(unit_t) min_dt
            close(unit_t)
         endif

         level_len = 0
         current => base%framework_next
         do while (current%lv .ne. 0)
            do i_level = 1, n_levels
               if (current%lv .eq. i_level) level_len(i_level) = level_len(i_level) + 1
            enddo
            current => current%framework_next
         enddo

         open(unit_datamap, file = trim(datamap_path) // 'data.t' // nt_char // '.cpu' // cpu_char // '.mapout',                   &
            form = 'unformatted')
         write(unit_datamap) n_blocks_vip_all
         write(unit_datamap) level_len
         write(unit_datamap) n_blocks(cpu_id_p1)

         current => base%next
         do while (current%lv .ne. 0)
            call block_output(current)
            current => current%next
         enddo

         close(unit_datamap)

         if (mod(nt, dn_save_state) .eq. 0 .and. cpu_id .eq. 0) then
            write(*, *)
            write(*, '(f6.2, a, i3, a, i3, a, f6.2, a)')                                                                           &
               t / t_end * 100d0, '% is completed, estimated remaining time:',                                                     &
               int(t_remain / 3600), 'h',                                                                                          &
               int(mod(t_remain, 3600.0d0) / 60), 'm',                                                                             &
               mod(t_remain, 60.0d0), 's'
            write(*, *)
         endif

         if (t .ge. ((n_save_state + 1) * dn_save_state - 1) * dt_output .and. save_state_flag .eq. 1) then
            n_save_state = n_save_state + 1
            call save_state(base)
            if (cpu_id .eq. 0) write(*,"(a, i6)") ' State save: ', n_save_state
         endif

         t_loop_pre = t_loop_pre * 100.0d0
         t_e = mpi_wtime ()
         time_output_one = time_output_one + t_e - t_s

         if (t .ge. t_end) exit

      endif

      max_et = 0.0d0

!----------------------------------------------------------------------------------------------------------------------------------|
!     AMR start
!----------------------------------------------------------------------------------------------------------------------------------|
      call amr (base)
!----------------------------------------------------------------------------------------------------------------------------------|
!     AMR stop
!----------------------------------------------------------------------------------------------------------------------------------|

      t = t + min_dt

      t_s = mpi_wtime ()

!      call check_info ('other')
      max_speed = 0.0d0
      max_te_ro = 0.0d0
      current => base%next
      do while (current%lv .ne. 0)
         call check_state (current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,         &
            current%x, current%y, current%z, nx_block, ny_block, nz_block)
         current => current%next
      enddo

      call mpi_allreduce(max_speed, max_speed_global, 1, mpi_real8, mpi_max, mpi_comm_world, error)
      max_speed = max_speed_global

      base%dt = cfl / (max_speed / base%dx + max_speed / base%dy + max_speed / base%dz)
      base%dt = min (cfl * minval((/ base%dx, base%dy, base%dz /)) ** 2 / (max_et + eps), base%dt)
      base%dt_hc = min (cfl * minval((/ base%dx, base%dy, base%dz /)) ** 2 / (thermal_amp * max_te_ro + eps), base%dt) * 1.0d0 + eps

      base%t = base%t + min_dt

      t_e = mpi_wtime()
      time_other_one = time_other_one + t_e - t_s

      ns = ns + 1
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  final step, free memory and output calculation time
!----------------------------------------------------------------------------------------------------------------------------------|
!   if (associated(base)) deallocate(base)
!   if (allocated(n_blocks)) deallocate(n_blocks)
!   if (allocated(level_len)) deallocate(level_len)
!   if (allocated(n_blocks_vip)) deallocate(n_blocks_vip)
!   if (allocated(blocks_vip_flag)) deallocate(blocks_vip_flag)
!   if (allocated(blocks_vip_flag_all)) deallocate(blocks_vip_flag_all)
!   if (allocated(hilbert_start)) deallocate(hilbert_start)
!   if (allocated(hilbert_end)) deallocate(hilbert_end)
!   if (allocated(i_gene)) deallocate(i_gene)
!   if (allocated(gene)) deallocate(gene)
!   if (allocated(n_blocks_update_send)) deallocate(n_blocks_update_send)
!   if (allocated(n_blocks_update_sendrecv)) deallocate(n_blocks_update_sendrecv)
!   if (allocated(blocks_update_info_send)) deallocate(blocks_update_info_send)
!   if (allocated(blocks_update_info_recv)) deallocate(blocks_update_info_recv)
!   if (allocated(n_neighs_reflux_sendrecv_x)) deallocate(n_neighs_reflux_sendrecv_x)
!   if (allocated(n_neighs_reflux_sendrecv_y)) deallocate(n_neighs_reflux_sendrecv_y)
!   if (allocated(n_neighs_reflux_sendrecv_z)) deallocate(n_neighs_reflux_sendrecv_z)
!   if (allocated(n_neighs_shareblock_sendrecv_x)) deallocate(n_neighs_shareblock_sendrecv_x)
!   if (allocated(n_neighs_shareblock_sendrecv_y)) deallocate(n_neighs_shareblock_sendrecv_y)
!   if (allocated(n_neighs_shareblock_sendrecv_z)) deallocate(n_neighs_shareblock_sendrecv_z)
!   if (allocated(n_neighs_exchange_sendrecv_x)) deallocate(n_neighs_exchange_sendrecv_x)
!   if (allocated(n_neighs_exchange_sendrecv_y)) deallocate(n_neighs_exchange_sendrecv_y)
!   if (allocated(n_neighs_exchange_sendrecv_z)) deallocate(n_neighs_exchange_sendrecv_z)
!   if (allocated(neighs_reflux_info_send_x)) deallocate(neighs_reflux_info_send_x)
!   if (allocated(neighs_reflux_info_send_y)) deallocate(neighs_reflux_info_send_y)
!   if (allocated(neighs_reflux_info_send_z)) deallocate(neighs_reflux_info_send_z)
!   if (allocated(neighs_reflux_info_recv_x)) deallocate(neighs_reflux_info_recv_x)
!   if (allocated(neighs_reflux_info_recv_y)) deallocate(neighs_reflux_info_recv_y)
!   if (allocated(neighs_reflux_info_recv_z)) deallocate(neighs_reflux_info_recv_z)
!   if (allocated(neighs_shareblock_info_send_x)) deallocate(neighs_shareblock_info_send_x)
!   if (allocated(neighs_shareblock_info_send_y)) deallocate(neighs_shareblock_info_send_y)
!   if (allocated(neighs_shareblock_info_send_z)) deallocate(neighs_shareblock_info_send_z)
!   if (allocated(neighs_shareblock_info_recv_x)) deallocate(neighs_shareblock_info_recv_x)
!   if (allocated(neighs_shareblock_info_recv_y)) deallocate(neighs_shareblock_info_recv_y)
!   if (allocated(neighs_shareblock_info_recv_z)) deallocate(neighs_shareblock_info_recv_z)
!   if (allocated(neighs_exchange_info_send_x)) deallocate(neighs_exchange_info_send_x)
!   if (allocated(neighs_exchange_info_send_y)) deallocate(neighs_exchange_info_send_y)
!   if (allocated(neighs_exchange_info_send_z)) deallocate(neighs_exchange_info_send_z)
!   if (allocated(neighs_exchange_info_recv_x)) deallocate(neighs_exchange_info_recv_x)
!   if (allocated(neighs_exchange_info_recv_y)) deallocate(neighs_exchange_info_recv_y)
!   if (allocated(neighs_exchange_info_recv_z)) deallocate(neighs_exchange_info_recv_z)
!   if (allocated(neighs_reflux_send_x)) deallocate(neighs_reflux_send_x)
!   if (allocated(neighs_reflux_send_y)) deallocate(neighs_reflux_send_y)
!   if (allocated(neighs_reflux_send_z)) deallocate(neighs_reflux_send_z)
!   if (allocated(neighs_reflux_recv_x)) deallocate(neighs_reflux_recv_x)
!   if (allocated(neighs_reflux_recv_y)) deallocate(neighs_reflux_recv_y)
!   if (allocated(neighs_reflux_recv_z)) deallocate(neighs_reflux_recv_z)
!   if (allocated(neighs_shareblock_send_x)) deallocate(neighs_shareblock_send_x)
!   if (allocated(neighs_shareblock_send_y)) deallocate(neighs_shareblock_send_y)
!   if (allocated(neighs_shareblock_send_z)) deallocate(neighs_shareblock_send_z)
!   if (allocated(neighs_shareblock_recv_x)) deallocate(neighs_shareblock_recv_x)
!   if (allocated(neighs_shareblock_recv_y)) deallocate(neighs_shareblock_recv_y)
!   if (allocated(neighs_shareblock_recv_z)) deallocate(neighs_shareblock_recv_z)
!   if (allocated(neighs_exchange_send_x)) deallocate(neighs_exchange_send_x)
!   if (allocated(neighs_exchange_send_y)) deallocate(neighs_exchange_send_y)
!   if (allocated(neighs_exchange_send_z)) deallocate(neighs_exchange_send_z)
!   if (allocated(neighs_exchange_recv_x)) deallocate(neighs_exchange_recv_x)
!   if (allocated(neighs_exchange_recv_y)) deallocate(neighs_exchange_recv_y)
!   if (allocated(neighs_exchange_recv_z)) deallocate(neighs_exchange_recv_z)
!   if (allocated(neighs_shareblock_send_x_hc)) deallocate(neighs_shareblock_send_x_hc)
!   if (allocated(neighs_shareblock_send_y_hc)) deallocate(neighs_shareblock_send_y_hc)
!   if (allocated(neighs_shareblock_send_z_hc)) deallocate(neighs_shareblock_send_z_hc)
!   if (allocated(neighs_shareblock_recv_x_hc)) deallocate(neighs_shareblock_recv_x_hc)
!   if (allocated(neighs_shareblock_recv_y_hc)) deallocate(neighs_shareblock_recv_y_hc)
!   if (allocated(neighs_shareblock_recv_z_hc)) deallocate(neighs_shareblock_recv_z_hc)
!   if (allocated(neighs_exchange_send_x_hc)) deallocate(neighs_exchange_send_x_hc)
!   if (allocated(neighs_exchange_send_y_hc)) deallocate(neighs_exchange_send_y_hc)
!   if (allocated(neighs_exchange_send_z_hc)) deallocate(neighs_exchange_send_z_hc)
!   if (allocated(neighs_exchange_recv_x_hc)) deallocate(neighs_exchange_recv_x_hc)
!   if (allocated(neighs_exchange_recv_y_hc)) deallocate(neighs_exchange_recv_y_hc)
!   if (allocated(neighs_exchange_recv_z_hc)) deallocate(neighs_exchange_recv_z_hc)
!   if (allocated(n_neighs_exchange_sendrecv_fix_y)) deallocate(n_neighs_exchange_sendrecv_fix_y)
!   if (allocated(n_neighs_exchange_sendrecv_fix_z)) deallocate(n_neighs_exchange_sendrecv_fix_z)
!   if (allocated(neighs_exchange_info_send_fix_y)) deallocate(neighs_exchange_info_send_fix_y)
!   if (allocated(neighs_exchange_info_send_fix_z)) deallocate(neighs_exchange_info_send_fix_z)
!   if (allocated(neighs_exchange_info_recv_fix_y)) deallocate(neighs_exchange_info_recv_fix_y)
!   if (allocated(neighs_exchange_info_recv_fix_z)) deallocate(neighs_exchange_info_recv_fix_z)
!   if (allocated(bnd_exterior)) deallocate(bnd_exterior)
!   if (allocated(bnd_interior)) deallocate(bnd_interior)
!   if (allocated(bnd_interior_pos)) deallocate(bnd_interior_pos)
!   if (allocated(thr_va0)) deallocate(thr_va0)
!   if (allocated(x_total)) deallocate(x_total)
!   if (allocated(ro_total_x)) deallocate(ro_total_x)
!   if (allocated(mx_total_x)) deallocate(mx_total_x)
!   if (allocated(my_total_x)) deallocate(my_total_x)
!   if (allocated(mz_total_x)) deallocate(mz_total_x)
!   if (allocated(bx_total_x)) deallocate(bx_total_x)
!   if (allocated(by_total_x)) deallocate(by_total_x)
!   if (allocated(bz_total_x)) deallocate(bz_total_x)
!   if (allocated(en_total_x)) deallocate(en_total_x)
!   if (allocated(y_total)) deallocate(y_total)
!   if (allocated(ro_total_y)) deallocate(ro_total_y)
!   if (allocated(mx_total_y)) deallocate(mx_total_y)
!   if (allocated(my_total_y)) deallocate(my_total_y)
!   if (allocated(mz_total_y)) deallocate(mz_total_y)
!   if (allocated(bx_total_y)) deallocate(bx_total_y)
!   if (allocated(by_total_y)) deallocate(by_total_y)
!   if (allocated(bz_total_y)) deallocate(bz_total_y)
!   if (allocated(en_total_y)) deallocate(en_total_y)
!   if (allocated(z_total)) deallocate(z_total)
!   if (allocated(ro_total_z)) deallocate(ro_total_z)
!   if (allocated(mx_total_z)) deallocate(mx_total_z)
!   if (allocated(my_total_z)) deallocate(my_total_z)
!   if (allocated(mz_total_z)) deallocate(mz_total_z)
!   if (allocated(bx_total_z)) deallocate(bx_total_z)
!   if (allocated(by_total_z)) deallocate(by_total_z)
!   if (allocated(bz_total_z)) deallocate(bz_total_z)
!   if (allocated(en_total_z)) deallocate(en_total_z)
!   if (allocated(unit_cell)) deallocate(unit_cell)

   t_stop = mpi_wtime()

   t_wall = t_stop - t_start

   if (cpu_id .eq. 0) then
      write(*, *)
      write(*, '(a)') '|-----------------|'
      write(*, '(a)') ' Program finished.'
      write(*, '(a)') '|-----------------|'
      write(*, *)
      write(*, '(a)') '|--------------------------------------------|'
      write(*, '(a, i3, a, i3, a, f6.2, a)') ' Total elapsed time:',                                                              &
         int(t_wall/3600.), 'h', int (mod (t_wall, 3600.) / 60), 'm', mod (t_wall, 60.), 's'
      write(*, '(a)') '|--------------------------------------------|'
      write(*, *)
      write(*, '(a)') '|--------------------------------------------------------------------------|'
      write(*, '(a)') ' Calculation costs:'
      write(*, '(a, f6.2, a)') ' Regrid:', time_regrid / t_wall * 100., '%'
      write(*, '(a, f6.2, a)') ' Integration:', time_integrate / t_wall * 100., '%'
      write(*, '(a, f6.2, a)') ' Data output:', time_output / t_wall * 100., '%'
      write(*, '(a, f6.2, a)') ' Boundary condition:', time_boundary / t_wall * 100., '%'
      write(*, '(a, f6.2, a)') ' Other (initial condition, data check, calculate minimum dt):', time_other / t_wall * 100., '%'
      write(*, '(a)') '|--------------------------------------------------------------------------|'
      write(*, *)
   endif

   if (cpu_id .eq. 0) then
      open(unit_par, file = trim(datamap_path) // 'params.mapout', position = 'append')
      write(unit_par, '(a24, t32, i3, a, i3, a, f6.2, a)') ' Total elapsed time:',                                                 &
         int(t_wall/3600.), 'h', int (mod (t_wall, 3600.) / 60), 'm', mod (t_wall, 60.), 's'
      write(unit_par, '(a24, t32)') ' Calculation costs:'
      write(unit_par, '(a24, t32, f6.2, a)') ' Regrid:', time_regrid / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Integration:', time_integrate / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Data output:', time_output / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Boundary condition:', time_boundary / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Other:', time_other / t_wall * 100., '%'
      close(unit_par)
   endif

   call mpi_finalize(error)

!----------------------------------------------------------------------------------------------------------------------------------|
   stop
end program main
